###############################################################################   
#### Replication Materials                                                 #### 
#### Kim, Nakka, Gopal, Desmrais, Mancinelli, Harden, Ko, Boehmke. 2021.   ####
#### Attention to the COVID-19 pandemic on Twitter:                        ####
#### Partisan differences among U.S. state legislators                     ####
#### Legislative Studies Quarterly                                         ####
###############################################################################  


###############################################################################
################################### Set Up ####################################
###############################################################################

# packages -------------------------

lapply(c('readr', 'stargazer', 'texreg', 'xtable', 'lme4', 'glmmTMB','plm',
         'lmtest', 'interplot','dummies','effects', 'ggthemes','caret', 
         'ggeffects', 'scales','ggrepel','sandwich','reshape2'), 
       require, 
       character.only = TRUE)

# read data set for regression -------------------------

path_data <- '/Users/taegyoon/Google Drive/spap_state/spap_state_attention/data/' 

df <- read_csv(paste0(path_data, "spap_state_attention_regression.csv")) 
df <- subset(df, select = -c(X1) ) 
df_dem <- df[which(df$party == 'D'), ] 
df_rep <- df[which(df$party == 'R'), ] 


# set the panel structure -------------------------

plm_df_dem <- pdata.frame(x = df_dem, 
                          index = c('user.screen_name', 'week'))
plm_df_rep <- pdata.frame(x = df_rep, 
                          index = c('user.screen_name', 'week'))


###############################################################################
#### Table S12: Regression with Dem-legislator Observations (legislator FE) ###
###############################################################################

# fit models -------------------------

legis_fe_state_dem <- plm(covid_relevant_1_log ~ 
                               state_case_pop + 
                               state_death_pop +
                               state_covid_policy + 
                               week_new +
                               week_new_2 +
                               week_new_3,
                             data = plm_df_dem,
                             index = c('user.screen_name', 'week'),
                             model = 'within',
                          effect = 'individual')

legis_fe_national_dem <- plm(covid_relevant_1_log ~ 
                                  national_case_pop + 
                                  national_death_pop +
                                  state_covid_policy + 
                                  week_new +
                                  week_new_2 +
                                  week_new_3,
                                data = plm_df_dem,
                                index = c('user.screen_name', 'week'),
                                model = 'within', 
                             effect = 'individual')

legis_fe_all_dem <- plm(covid_relevant_1_log ~ 
                             state_case_pop + 
                             state_death_pop +
                             state_covid_policy + 
                             national_case_pop + 
                             national_death_pop +
                             week_new +
                             week_new_2 +
                             week_new_3,   
                           data = plm_df_dem,
                           index = c('user.screen_name', 'week'),
                           model = 'within', 
                        effect = 'individual')


# adjust SE for serial correlation -------------------------

legis_fe_state_dem_adj <- coeftest(legis_fe_state_dem, 
                                   vcovHC(legis_fe_state_dem, 
                                          method='arellano')) 
legis_fe_national_dem_adj <- coeftest(legis_fe_national_dem, 
                                      vcovHC(legis_fe_national_dem, 
                                             method='arellano'))
legis_fe_all_dem_adj <- coeftest(legis_fe_all_dem, 
                                 vcovHC(legis_fe_all_dem, 
                                        method='arellano'))


###############################################################################
######################## Generate Latex Code : Table S13 ######################
###############################################################################

# get information for the last three rows -------------------------

texreg(
  list(
    legis_fe_state_dem,
    legis_fe_national_dem,
    legis_fe_all_dem),
  digits = 3,
  caption = 'Panel Regression Model (Legislator Fixed Effect): 
  Democratic Legislators, Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)'
                        )
  )

# get information for betas and corresponding standard errors -------------------------

texreg(
  list(
    legis_fe_state_dem_adj,
    legis_fe_national_dem_adj,
    legis_fe_all_dem_adj),
  digits = 3,
  caption = 'Panel Regression Model (Legislator Fixed Effect): 
  Democratic Legislators, Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)'
                        )
  )


###############################################################################
#### Table S12: Regression with rep-legislator Observations (legislator FE) ###
###############################################################################

# fit models -------------------------

legis_fe_state_rep <- plm(covid_relevant_1_log ~ 
                            state_case_pop + 
                            state_death_pop +
                            state_covid_policy + 
                            week_new +
                            week_new_2 +
                            week_new_3,
                          data = plm_df_rep,
                          index = c('user.screen_name', 'week'),
                          model = 'within',
                          effect = 'individual')

legis_fe_national_rep <- plm(covid_relevant_1_log ~ 
                               national_case_pop + 
                               national_death_pop +
                               state_covid_policy + 
                               week_new +
                               week_new_2 +
                               week_new_3,
                             data = plm_df_rep,
                             index = c('user.screen_name', 'week'),
                             model = 'within', 
                             effect = 'individual')

legis_fe_all_rep <- plm(covid_relevant_1_log ~ 
                          state_case_pop + 
                          state_death_pop +
                          state_covid_policy + 
                          national_case_pop + 
                          national_death_pop +
                          week_new +
                          week_new_2 +
                          week_new_3,   
                        data = plm_df_rep,
                        index = c('user.screen_name', 'week'),
                        model = 'within', 
                        effect = 'individual')


# adjust SE for serial correlation -------------------------

legis_fe_state_rep_adj <- coeftest(legis_fe_state_rep, 
                                   vcovHC(legis_fe_state_rep, 
                                          method='arellano')) 
legis_fe_national_rep_adj <- coeftest(legis_fe_national_rep, 
                                      vcovHC(legis_fe_national_rep, 
                                             method='arellano'))
legis_fe_all_rep_adj <- coeftest(legis_fe_all_rep, 
                                 vcovHC(legis_fe_all_rep, 
                                        method='arellano'))


###############################################################################
######################## Generate Latex Code : Table S13 ######################
###############################################################################

# get information for the last three rows -------------------------

texreg(
  list(
    legis_fe_state_rep,
    legis_fe_national_rep,
    legis_fe_all_rep),
  digits = 3,
  caption = 'Panel Regression Model (Legislator Fixed Effect): 
  repocratic Legislators, Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)'
                        )
  )

# get information for betas and corresponding standard errors -------------------------

texreg(
  list(
    legis_fe_state_rep_adj,
    legis_fe_national_rep_adj,
    legis_fe_all_rep_adj),
  digits = 3,
  caption = 'Panel Regression Model (Legislator Fixed Effect): 
  repocratic Legislators, Population Normalized + Three-fold Party Variable + Week Polynomial + Logged DV (1 added)',
  custom.coef.names = c('State New Cases (per 10k)',
                        'State New Deaths (per 10k)',
                        'State COVID-19 Policies',
                        'Week',
                        'Week (quadratic)',
                        'Week (cubic)',
                        'National New Cases (per 10k)',
                        'National New Deaths (per 10k)'
                        )
  )